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CONVERGENCE OF DIRECT RECURSIVE ALGORITHM FOR 
IDENTIFICATION OF PREISACH HYSTERESIS MODEL WITH 

STOCHASTIC INPUT 

D. RACHINSKII* AND M. RUDERMANt 

Abstract. We consider a recursive iterative algorithm for identification of parameters of the 
Preisach model, one of the most commonly used models of hysteretic input-output relationships. 
The classical identification algorithm due to Mayergoyz defines explicitly a series of test inputs that 
allow one to find parameters of the Preisach model with any desired precision provided that (a) such 
input time series can be implemented and applied; and, (b) the corresponding output data can be 
accurately measured and recorded. Recursive iterative identification schemes suitable for a number 
of engineering applications have been recently proposed as an alternative to the classical algorithm. 
These recursive schemes do not use any input design but rather rely on an input-output data stream 
resulting from random fluctuations of the input variable. Furthermore, only recent values of the 
input-output data streams are available for the scheme at any time instant. In this work, we prove 
exponential convergence of such algorithms, estimate explicitly the convergence rate, and explore 
which properties of the stochastic input and the algorithm affect the guaranteed convergence rate. 

Key words. Identification problem, recursive algorithm, exponential convergence rate, model 
of hysteresis, input-output operator, stochastic input. 
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1. Introduction. Preisach model [22], with its roots originated in magnetism 
[18,19], is well-suitable for describing various rate-independent hysteresis phenomena 
in a closed input-output analytic form. The Preisach hysteresis operator y(k) = 
( T-L[x\)(k ) models an inpiut-output relationship with an erasable memory of previous 
input states x{j ), where j < k. The input-output time series ( x,y)(k ) constitute 
monotonic piecewise continuous hysteresis curves, see Figure 1 (right). 



Fig. 1. Preisach density function p(a,/3) over the plane (ex, (3) (left) and input-output curves 
y = ('H[x])(t) of the Preisach operator (right). 
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According to the representation theorem of Mayergoyz, the Preisach hysteresis 
model can be fitted to any set of input-output data that possesses two properties: 
the so-called wiping-out of local input extremum values and congruency of any two 
hysteresis loops created by the same two consecutive input extrema; see [18] for more 
details of the formalism and properties of the Preisach hysteresis model. 

The use of models of hysteresis in magnetism, as well as in other natural, engineer¬ 
ing, and social sciences, confronts with problems of identification of parameters, which 
requires to obtain a sufficient number of observations of hysteresis system states [2,8]. 
The phenomenological nature of the Preisach model, which is based on superposition 
of multiple elementary hysteresis operators called non-ideal relays (switches), and the 
representation theorem of Mayergoyz allow to apply the model to various hysteretic 
systems using a phenomenological argument or the black box approach without com¬ 
plex analysis of such systems based on first principles [3,7,10,11,15-17,20]. Impor¬ 
tantly, the representation theorem provides the Preisach model with straightforward 
means for identification of parameters from dedicated input-output data. Here it 
should be recalled that the Preisach model is entirely parameterized by the so-called 
Preisach density function p = p(a, 0) defined over the {a,0) plane, as schematically 
shown in Figure 1 (left). The a and 0 coordinates, with 0 > a, are the ‘on’ and 
‘off’ switching thresholds of a single relay operator, and the p-function represents the 
weighted ( a , /^-distribution of these relays over the possible threshold values. Each 
relay switches from state 0 (‘off’) to state 1 (‘on’) and backwards in response to 
variations of the same input time series called the input x(k) of the Preisach model; 
the output time series y(k ) of the model is defined as the sum of states of all the 
relays. The identification algorithm of Mayergoyz requires to discretize the Preisach 
model and assumes a sufficient number of measured first-order descending (FOD) 
input-output curves [18]. The higher the required accuracy of hysteresis modeling is, 
the smaller the discrete step on the (a,/3) plane should be (cf. with Figure 1, left), 
and the more accurate FOD curve measurements are required. The volume of data, 
which should be first recorded before the direct Mayergoyz identification method can 
be applied for calculation of the density function p, rapidly increases with increasing 
target accuracy. 

From this point of view, the motivation for this work is to explore the possibili¬ 
ties of applying recursive/iterative identification methods for which only the previous 
(or, a few recent) values of the input-output data stream should be available at any 
instant. Recursive (also known as on-line) identification methods can be required 
when dealing with time-varying processes, i.e., where the hysteresis behavior changes 
(usually slowly) with time depending on certain non-controllable operational or en¬ 
vironmental conditions, such as varying ambient temperature, force fields, temporal 
relaxation, etc. The applicability of the direct Mayergoyz identification method can 
also be limited by a strong uncontrollable noise capable of corrupting the designed 
test input signal. Furthermore, in various engineering applications which use hys¬ 
teresis models, e.g., for process monitoring or control, a simplified identification with 
no special (expert) knowledge as for design of dedicated input processes (such as a 
process creating the FOD curves) can be required. Here an algorithm which is capa¬ 
ble of using data generated by a stochastic input processes (with certain well-defined 
properties) may be advantageous. 

Surprisingly, little number of recursive methods suitable for identification of pa¬ 
rameters of the Preisach hysteresis model has been analyzed theoretically or experi¬ 
mentally. A gradient-based recursive identification scheme has been proposed in [28] 
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and experimentally evaluated on a magnetostrictive actuator system which exhibits 
hysteresis between the relative displacement and the actuator current. The error 
convergence for this scheme has been estimated for relatively rough discretization 
of the (ct,/3) plane (up to 25 steps in the input domain). Another direct recursive 
identification method has been proposed in [24], This method updates the values of 
the Preisach density function in the so-called switching region of the (a, 0) plane, 
which defines an output increment in response to a random input increment at each 
iteration step. This method has been tested and evaluated numerically in the con¬ 
text of modeling the hysteretic constitutive relationship between the magnetic field 
and magnetization in a ferromagnetic material. However no rigorous analysis of the 
convergence and its properties, such as the convergence rate, has been done insofar. 

In this paper, we prove the exponential convergence of the recursive identification 
algorithm proposed in [24], We first show that the error of the algorithm monoton- 
ically decreases (Section 2) and establish an estimate for the rate of convergence for 
any input time series (Section 3). Our main result (Section 4) establishes that the 
error Ek almost surely converges to zero exponentially, ||EA,|| < Ce~ A:A ||i?o||, with the 
rate satisfying an explicit estimate A > y(No, (T ), (1 /T)) > 0 if the input time series 
generates a regular Markov chain in the space of states of the Preisach model (the 
latter condition is generic and is satisfied for a wide class of input Markov processes). 
Here Nq is a function of the number of nodes in the domain of the density function on 
the (a, 0) plane, T is a certain stopping time associated with the input process, and 
(T), (1 /T) are the mean values of T and its inverse. Essentially, T is the length of the 
time interval during which fluctuations of the input randomly create all the segments 
of all the FOD input-output curves prescribed by the classical Mayergoyz identifica¬ 
tion algorithm. If all this information was recorded, a complete identification of the 
density function was achievable from the input-output data collected within T time 
steps. The recursive algorithm considered here is not capable of recovering the den¬ 
sity function exactly as it updates the approximation to the density function at each 
time step without recording past input-output data. However, we guarantee that the 
error decreases by a certain percentage over T time steps, UE/c+tII < e - c /( iV oT)||_g’ fe ||. 
It is worth noting that the random time T is a characteristic of the input process. 
Overall, our estimates suggest that the smaller the mean value of T is and the less 
nodes the grid has (smaller N 0 ), the faster is the guaranteed exponential convergence 
of the recursive algorithm. We test these conclusions by considering two examples of 
random input processes and comparing their rates of convergence in Section 5. 

2. Monotonic decrease of the error. Let p = p(p) be a density function of 
the Preisach operator where p = (a , 0) denotes a point in the half-plane n = {p = 
( a,/3 ) : a < /3}. Each point p represents a non-ideal relay with the threshold values 
a and 0 Denote by Sk = S(tk) the subset of the half-plane n where relays are ‘on’ 
(in state 1) at a moment tk = k > 0 (we normalize the time step to 1); the other 
relays are ‘off’ (in state 0). Then the output yk = y{tk) of the Preisach model at this 
moment equals 



where dp = dad0 The set Sk is called the state of the Preisach model. The state 
Sk and the output yk change in response to variations of the input Xk = x(tk ) of the 
model as the non-ideal relays switch between states 0 and 1. For discrete time input 
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series, which is the case considered here, 


S k U {(a, /?) € II : /3 < x k+ i} if x k +i > x k 
S k +i = { 4\{(a,(3)eII:«>i H i} if x k +i < x k 
S k if x k -\-i — x k 


This rule defines the time series of the state and output, S k and y kl for any given input 
sequence x k , k > 0 and a given initial state So (see, for example, [4-6,9,12-14,21,23,29] 
for more details). 

In the following identification algorithm we use an estimate p k = p k (p ), p € II of 
the density function, which is updated at each moment t k . The estimate y k = y{t k ) 
of the output and the error e k = y k — y k between the observed and modeled output 
values equal, respectively, 


Ilk 






( P(P ) - Pk(p)) dp. 


Comparing states S k and S k+ 1, denote f = S k +i \S k . This is the set of those relays 
that switch on during the time step from t k to t k + 1. Similarly, consider the set 
= S k \ S k + i of relays that switch off during the same time step. Then, over this 
time step, the error between the observed output and the output modeled with the 
density function p k receives the increment 


Ae k 


/ {pip) ~ Pk{p )) dp- (p{p) - p k (p)) dp. 

Jn+ Jn k 


Following [24], we update the estimate of the true density function p according 
to the formula 


(2.1) p k+1 (p) = p k (p) + r k {p) with r k (p) = < 


Ae fc 

Asi. 


if p flfc = U 

if p € Ilfc 


Ae fc if „ p O- 
A n 11 


where Ao, k denotes the area of the switching region Q k . 

Let us first show that the mean square distance between the exact and estimated 
density is monotonically non-increasing with iterations. 

Lemma 2.1. For any initial estimate po, the mean square distance 


II Pk - p\\ = ~ PiP )) 2 d P^j 


between p k and p is a non-increasing function of k. 
Proof. According to the update rule, 


\\p k +i-p\\ 2 = / {pkip) + r k (p) - p{p)) 2 dp 
J n 


iPk{p) - p{p)) 2 dp + 


r t ( P ) dp + 2 r k (p) ( p k {p) - p{p)) dp 
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= \\p k ~p\\ 2 + 


(Ae fc ) 2 


2Ae fc 

A n k 


/ ihip) - Pi?)) dp 



pip)) dp 


= \\Pk ~ p \\ 2 + 


(Aefc ) 2 


2(Ae fc ) 2 


llPfc — pH 2 


(Ae fc ) 2 


□ 


We see that 

(2-2) ||pfc +1 -p|| 2 = ||pfc-p|| 2 -^^, 

hence the mean square distance between the exact and estimated Preisach density 
functions not only never increases, but actually strictly decreases each time the output 
increment error Ae is non-zero. The aim of the following argument is to show that, 
for reasonable inputs, the error Aefc is non-zero “often enough” (as long as p k ^ p) 
to ensure the exponential convergence p k —> p. 

3. Estimation of the rate of convergence. From now on, we will consider the 
discretized version of the Preisach model. We will assume that the input is a discrete 
Markov process X k with N + 1 states Xi = is, where e = 1/N and i = 0 ,N. 
In this discrete setting, the set S of all possible states S of the model is also finite. 
Furthermore, we can assume without loss of generality that the density p (as well as 
its approximations p k ) is concentrated at the finite set of nodes 


(3.1) {(a, 0) = {{t - l/2)e, (m - l/2)e) : 1 < t < m < N} 


that form a uniform rectangular mesh within the right triangle 

A supp = {(a, 0) ■ 0 < a < /3 < 1} 

of the half-plane II; these nodes will be denoted pi = (a,,^*), where i € I with the 
total number of nodes #(/) = M = N(N + l)/2. In other words, we consider the 
density functions of the form p(p) = e 2 )T) i=1 p{pi)5{p — pi). It is convenient to rescale 
the mesh step to unity by setting 

Vk = Y Y PkiPi), Aefc = Y iPiPi)~Pk{Pi))~ Y (P(Pi)~Pk(Pi)) 

i€i iei iei iei 

Pi^Sk Pi^Sk Pi^Q^ 


and, simultaneously, setting 

An k = N k := #(* £ I : Pi € fife) 

in the update rule (2.1) with N k denoting the number of the mesh nodes in the 
switching region fife. We note that for this discrete model, at each time step either 
fl+ = 0, flfc = fl)7 or fl)7 = 0, fife = f . With this notation, formula (2.2) can be 
written as 

(3.2) ||£fc+i|| 2 = H^fcll 2 - N k {E k+l {p) - Efc(p)) 2 , Vp e fife, 

where E k is the error 


E k {p) = pkip) - Pip) 
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1 /9 

and H^fell = (E ie j Ek(,Pi) 2 ) ■ Relationship (3.2) can be also written as 

(3-3) H^fc+ill 2 = ||£fc|| 2 - \\E k ~ E k+1 \\ 2 . 

We define the rate of convergence of the iterative scheme (2.1) as 


(3.4) 


A = — lim sup — In 

fc—»oo & 


\\E k \\ 

ll^oll 


1 

lim sup — y In 

k—y oo ^ 1 

i—l 


m\ 


In order to estimate the rate of convergence of the identification algorithm, we 
will use a certain stopping time T associated with the input process. With a given 
realization of the input process X k we associate successive time intervals 


1 < fc < Ti, T\ + 1 < k < Ti + T 2 , 


T 1 +T 2 + l<k<T 1 + T 2 + T 3 ,... 


During each time interval we require the process to create a set of Preisach states 
satisfying the following properties. We require that for every node pi = (ai,fa) there 
would be two states S(t k i) and S(t k > +1) achieved at two successive moments 1 within 
the time interval pi + ■ ■ ■ + T k _\ + 1, Tf + • • ■ + T k \ such that 

(5(t fc 0 \ S(t k . + 1)) U ( S(t k , + 1) \ S(t k ,)) = Aj 


where 

(3.5) A i = {(a, fa) : a, - e/2 < a < /3 < ft + e/2}. 


As soon as this requirement is satisfied for all the nodes (), i £ I, of the mesh, 
the new time interval [Ti + ■ ■ ■ + T k + 1, Ti + • —K Tfe+i] starts. The set (3.5) is a right 
triangle with the vertex (a, — e/2, + e/2) and the hypothenuse on the line a = /3. 

Therefore, we will simply say that the input creates all the possible triangles on the 
half-plane II within each time interval pi + • • • + T k _ 1 + 1, Tf + ■ • ■ + T k \. 

For example, this requirement will be satisfied if for every ( ai,/3i ) there will be 
three successive time moments t k = k , t k +i = k + 1 and t k + 2 = k + 2 during the time 
interval pi + • ■ • + T k + 1, Ti H-b Tfe+i] such that either 

(3.6) x{t k ) = /3j + e/2, x(t k+1 ) = - e/2, x(t k+2 ) = fa + e/2 

for some fa > fa, or 

(3.7) x(t k ) = aj - e/2, x{t k+1 ) = fa + e/2, x{t k+2 ) = a* - e/2 

for some aj < Oj. A less restrictive requirement is that the same holds but the 
moments t k , t k +\ are not necessarily successive and the input stays between — e/2 
and fa + e/2 for t k < t < t k+ \ in the first case (between aj — e/2 and fa + e/2 for 
t k <t < t k+ i in the second case). 

Theorem 3.1. The rate of convergence (3.4) of the iterative scheme (2.1) satis¬ 
fies the estimate 


(3.8) 


A > 


m( k ) 

-lim inf — 

32JVo fc—foo k 

i= 1 


1 

1 + Ti 


1 It is easy to extend the algorithm and results to the case when these moments are not necessarily 
successive. We consider successive moments for simplicity. 
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where m = m(k) denotes the largest integer such that k > T\ + ■ ■ ■ + T m and 


(3.9) 


N 0 


N(N + l) 2 i N 2 (N 2 — 1) 

4 + 24 


In the next section, we will show that the limit in (3.8) is positive almost surely 
and can be found from the stationary distribution of the input process under natural 
assumptions. 

Proof of Theorem 3.1. Define tq = 0, r m = T\ + • • • + T m and assume that the 
input creates all the possible triangles during each time interval [rfc_i,Tfc]. Set 

^=Pr m _ 1 || 2 -||^ m || 2 . 


According to formula (3.3), 

Tm~ 1 

*= E \\Ek+i-E k \\ 2 . 

k—Tm — 1 

We split this sum into M sums, where M is the total number of nodes: 

M Tm~ 1 

= E \\E k+ i-E k \\ 2 

2=1 /c—T m _ i 


where the coefficients 9i satisfy 


0i + • • ■ + 0m — I- 


Now, we replace the upper limit r m — 1 in the internal sum by the moment U £ 
[r m _i,r m — 1] when the switching area during the time step from ti to U + 1 is the 
triangle Aj = {(a,/3) : at — e/2 < a < /3 < + e/2}. In this way, we obtain the 

estimate 

M ti 

^E* E ||^fc+i — Ek \\ 2 ■ 

2=1 k=T m —1 


At the time step from ti to ti + 1 we have 

\\E ti+1 ~ E ti \\ 2 = 

where A\ = Et(p) denotes the sum of E t (p) over the triangle A i at the moment 

t and Ni is the number of nodes in the triangle. Hence, 


a > 


M 




{A'i? 

Nt 


M ti -1 

+E^ E ii^+i-^ii 2 - 

2=1 k—Tm —1 


Using the inequality 

E \\E k+ i-E k \\ 2 >—± 


k —r m — i 


ti 7~m—l 


( Ek+i — Ek) 


k=Tm — 1 
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we obtain 




k—Tm — 1 


and therefore 


M 


^E< 


7m-l 


A4‘) 2 , 11^. 
V JVi 


T„ 


T„ 


Let us consider one term in this sum: 
(3.10) 

(A ‘*) 2 , || E ti -E Tm _ x f_ 

' N, 


Ni 


Trr 


> ^ (E ^(p)) + ^- E (^(p) - ^-i(p )) 2 


t^peAi j P eAi 

The right hand side of this inequality can be rewritten as 


1 

Ni 


E (-^T-m-l (p) + Xp) I + TFT E 


^.peAi 


/yi Xp 

m p6Ai 


where \p = E u(p ) — E T m ~\ (p). Considering the minimization problem 


• • ■ >Cjv») = TT ( E C^Tm-i( p) + £p) ) + E 5 


ypeAi 


> ^P 
m pGAi 


with respect to the variables £ p , we obtain the equations 


E ( E r m - 1 (p) + Cp) ] + Jr~ — 0) P S Aj, 




ipeAj 


for the extremum, hence the minimum is achieved at the point 

T 

-L m. 


6 = ■ • ■ = (n, = - 


Ni(l+T m ) 


E E rm-AP) 


pe Ai 


and the global minimum value equals 

Vmin — 


Aj(l + T m ) 


E E r m -AP) I = 


(A 


T m- 1\2 


kpGAi 


iVi(i + r m )' 


Thus, this expression is the lower estimate for the the right hand side of inequality 
(3.10) and we conclude that 
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We rearrange the sum in the right hand side as follows: 

. i 0Mrr | wi- 1 ) 2 , e^ry \ 

- 4(1 + T m ) Ni Nj + N k N t j 

where j = j(i), k = k(i), £ = £{i) are defined so that A j = {(a, P) : on+e/2 < a < p < 
Pi + e/2}, A fe = {(a, f3) : a, - e/2 < a < ft < Pi - e/2}, At = {(a, P) : a, + e/2 < a < 
P < pi—e/2} (for the nodes (a,, Pi) on the diagonal P = a we set Aj = A k = At = 0). 
We choose 



with 

M N N 

N 0 = j2 N i = pJ2 k ( k + X )( N ~ k + l ) = 2 + x ) k +- fc3 ) 

i= 1 fc=l /c=l 


7V(iV + l ) 2 ( N 2 (N + 1)(27V + 1) iV 2 (iV + l ) 2 iV(7V + l ) 2 ( 7V 2 (7V 2 - 1) 
" 4 + 12 8 ~ 4 + 24 

where the total number of nodes is M = N(N + l)/2. With this choice of 6i , 

M 


a > 


1 


4(1 + T m )No 
We now note that for each node i 


£ ((A—) 2 + (A- 1 ) 2 + M-‘f + (4~‘> 2 ) ■ 


A T ~-' - A - Al" 1 + AJ "- 1 = E Tm _,(i). 

Minimizing the function W(£i, £ 2 , £ 3 , £ 4 ) = £ 2 + £ 2 + £3 + Cl under the constraint 
£1 + £2 + £3 + £4 = E Tm _ 1 (i), we see that min W = {E Tm _ 1 (i)) 2 / 4, hence 


M 


^=ll^ m _ 1 || 2 -||^ Tm || 2 > 


4(1 + T m )N 0 ^ 4 


^o^eo) 2 ii^.j 2 


16(1 + T m )No ' 


That is, 


ll^ m || 2 < 1 - 


1 


16(1 + T m )N 0 
As ln(l — x) < —x for 0 < x < 1, this implies 


p Tra _J 2 . 


_ln ll- E T m || +ln||-Er m _ 1 || > 


32(1 + T m )N 0 ' 


Summing, we obtain 


1 m 1 

-ln||S fc ||+hi||^||>—-X! 


32Nq . 1 + Ti 

i=i 


where we use the largest m = m(/c) such that k > T\ + • • • + T m . Therefore, the rate 
of convergence (3.4) of the iterative scheme satisfies (3.8). □ 
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4. Main result. As the input increment Xk+\ — Xf and the state Sk of the 
Preisach model at a moment k define the state Sk+i at the next moment, the input 
Markov chain Xk defines the Markov chain Sk in the state space S of the Preisach 
model describing the random evolution of the state. 

Now, let us consider two more auxiliary Markov chains defined by the input 

process. One is Yk = Sth _ \-T k - That is, any realization of the Markov chain Sk 

generates a realization of the auxiliary process Yk = S Tk where all the triangles (3.5) 

are created between any two moments Tfc-i = T\ H- \-Tk-i and Tk = Tl + ■ —b Tk 2 . 

The finite state space of the Markov chain Yk is the space S of states of the Preisach 
model. Here and henceforth we assume that the Markov chain Sk almost surely 
creates all the triangles (3.5) in finite time, hence the Markov chain Yk = S Tk is well 
defined. 

We will assume that properties of the input Markov chain Xk ensure that the 
Markov chain Yk = S Tk is regular (ergodic), that is some power of the transition 
probability matrix of the process Yk is strictly positive. 

Secondly, we consider an auxiliary Markov chain Zk = (Sth_ T k ,Tk) with a 

countable number of states. A pair ( S, T) belongs to the state space Z of this Markov 
chain if and only if there is a Preisach state S such that P(St = S, Ti = T|So = 
S) > 0. Let us observe that the transition probability P(Zk+ i = Z'\Zk = Z) from a 
state Z = (S, t) £ Z to a state Z' = (S' ,t') £ Z of the Markov chain Zk is, by its 
definition, independent of the component r of the initial state: 

(4.1) P(S Tl +-+T k+1 = S', T k+ i = t'| S Tl +-+T k = S, T k = t ) =: ps,S',r'- 

As Yk = S Tk is a regular Markov chain, for every Preisach state S £ S there is at 
least one positive integer T such that (S, T) £ Z. 

Lemma 4.1. Assume that the Markov chain Yk = Sth_ \-T k is regular. Then, 

the Markov chain Zk = (St h_ \-r k ,Tk) is (a) irreducible; (b) positive recurrent; and, 

(c) aperiodic. 

Proof. 

Irreducibility. Consider any pair of states Z = (S,T),Z = (S,T) £ Z. By 
definition of the state space Z, there is at least one state Z = (S,T) £ Z such that 
P(S’t 1 = S, T\ = T’lS'o = S) > 0. On the other hand, as the process Yk = S Tk 

is regular, there is a positive integer ko such that P(St 1 a _ hT ko = S|So = S) > 0. 

Combining these two estimates and using the time invariance of the Markov chain Sk, 
we obtain T > (St 1 +...+t ) co+1 = S,Tk 0 + 1 = T|So = S) > 0. But this is the probability 
that the process Zk reaches the state Z from the state Z in ko + 1 steps. As the 
probability is positive, we conclude that Zk is irreducible. 

Positive recurrence. For any state Z = (S, T) £ Z consider the state 
Z = (S,T) £ Z such that p := P(S't 1 = S, r I\ = T’lS'o = S) > 0. Consider again a 
sufficiently large integer k 0 such that the regularity of the process Yk = Sr k implies 
Pm '■= m i n {-F > (Sri+...+T fco = S|So = S) : S £ 5} > 0. We see that the probability 
that the process Zk starting from the state Z returns to this state after ko + 1 time 
steps, P(S Tl +-+T ko+ 1 = S, T feo+ i = T|S 0 = S), is greater or equal than p m p > 0, 

hence the Markov chain Zk is recurrent. Furthermore, P(St 1 h_ \-T ko+1 = S,Tk 0 +1 = 

T|So = S) > PmP for any S £ S. Therefore, the probability that the first return time 

2 According to the definition of the times T k , each moment r k is minimal in the sense that at 
least one triangle is not created between the moments r^_ i and — 1. 
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to the state Z, exceeds n(k 0 + 1) is less or equal than (1 — p m p) n for any integer n. 
Hence, the expected value of this first return time does not exceed the value 

fco ~I - 1 2fco+2 3/co+3 

^2 k + (l - p m p) ^2 k+(l-p m p) 2 ^2 k-\ - 

k—1 k=ko~\~2 /c—2 /cq+3 

As this sum is finite, the expected value of the first return time is finite, hence the 
process Zk is positive recurrent. 

Aperiodicity. Using the same notation as above, p = P(S = 5,4 = 
T\Sq = 5) > 0. From the regularity of the process Yfc = S Tk it follows that p{k) := 
P(S Tl + - +T k = ^l-^o = S) > 0 for any k > k 0 . Therefore, P(S Tl +-+T k+1 = S,T i = 
T|5o = 5) > 0, that is the process Zk returns to the state Z = (S,T) with positive 
probability after k time steps for any k > k 0 + 1. Hence the state Z is aperiodic. 
Since this is true for any Z € Z, the Markov chain Zk is aperiodic. □ 

Lemma 4.1 ensures that the finite Markov chain Yk has a strictly positive and 
unique stationary distribution 7ry. Furthermore, the Markov chain Zk also has a 
unique stationary distribution 7 tz and the ergodicity property 

^ 2=1 

holds for any initial state Z\ and any function / with a finite expected value 


E nz[f( Z )] < °0 

with respect to the distribution nz (see, for example, [27]). In particular, using 
f{Z) = T where Z = (5,T), we obtain 

1 k 

(4.2) -^T fc 4i^ z [T fc ]=:T* 


(the condition T* < oo will be established in the proof of Theorem 4.2 below). That 
is, the time average of the time intervals Tk converges to the mean T* of Tk with 
respect to the stationary distribution of the chain Zk ■ Similarly, 


(4.3) 


E 


k 4^ 1 + Tfc 
2=1 


E n 


1 


1 + 4 


=: T** 


Let us remark that the mean value T* is defined by 

T*= ^2 T ' ' ( 7r ^)s , ,T' 

(, S',r')ez 

where ( 7 tz)s',t' is the probability to find the process Zk in the state (Sk, Tk) = ( S ', r') 
according to the stationary distribution of this process. According to (4.1), for the 
stationary distribution, 

(tz)s',t' = ^2 ( 7r - z )<S',' r ' PS,S',r ' = E ns 'PS,S',t> 

(S,r)eZ sgs 


(4.4) 
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where Ilg = '^2 t {'^z)s,t and we sum over all the pairs (S', r) £ Z with the fixed S. 
Summing (4.4) with respect to r' over all pairs (S',r') £ Z with the fixed S', we 
obtain 


IIs' = ^2 P S,S'^S 
ses 

where Ps,s ' = PS,S' ,t‘ ■ As the stationary distribution ^ry of the Markov chain 

Yfe = Xt x -i_| -T k satisfies the same equation with the same transition probabilities 

Ps,S‘ and the same normalization condition, the uniqueness implies II = 7Ty. Hence, 
(4.4) is equivalent to 


( kz)s',t' = J^(t ty)s ■PS,S',r' 
ses 

and, furthermore, the mean T* = E nz [Tk] of the time Tk satisfies 

(4.5) E„m = Yl (ttf)s Y t '-PS,S',t'- 

ses (S',t')c2 


That is, we average the time required to create all the triangles (3.5) (that is, the 
transition time from Yj, to Yk+ 1 ) over the stationary distribution 7ry of the initial 
state Yk = S (regardless of the destination state Yk+i = S'). Similarly, 


(4.6) 


Ett 


1 


1 + T k 


ses 


PS,S', T 


1 + r' 

(, S',r')ez 


Now, we are ready to formulate the main result. 

Theorem 4.2. Assume that the Markov chain Yk = SViH _ \-T k induced by the 

input process Xk in the Preisach state space S is regular. Then the error of the itera¬ 
tive scheme (2.1) exponentially decreases and the rate of the exponential convergence 
(3.4) satisfies the estimate 


(4.7) 


A > 


Ett z 


1 

1 +Tk 


32N 0 -E vz [T k ] 


> 0 


where Nq is given by (3.9) and E nz [Tk], E^ 


1+T fc 


are defined by (4.5), (4.6). 


According to (4.7), the lower bound for the rate of the exponential convergence is 
controlled by the average time that the input process takes to create all the triangles 
on the half-plane n. One can assume that the shorter is this time, the faster is the 
convergence. We will test this conjecture with numerical examples in the next section. 
Also, the lower bound decreases as 0(N~ 4 ) with the increasing density of the mesh 
on the half-plane n. 


Proof of Theorem f.2. 


From Tk > 1 it follows that E^ z 


l 

1 +T k 


< oo, hence 


Lemma 4.1 implies relation (4.3) [27]. 

Let us show that also E^ z [Tk] < oo. Indeed, we assumed that the Markov chain 
Sk starting from any initial state So almost surely creates all the triangles (3.5) in 
finite time. Therefore, there is an integer r* such that the process Sk creates all the 
triangles (3.5) within the first r* time steps with a positive probability ps 0 > p* > 0 
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for any S 0 € S. Hence, the probability (4.1) satisfies ps,s',t' < (1 — P*) n for every 
t 1 > tit * and every S, S' £ S and consequently formula (4.5) implies E nz [T^] < oo. 

As the mean E nz [T*,] of the times T\ is finite, using Lemma 4.1 we conclude that 
(4.2) is valid. 

Finally, let us see that 


m(fc) 

(4.8) - ^ Ti 4 1 

i— 1 

where m = m(k) denotes the largest integer such that k > T) + • • • + T m . For any 
e > 0, the estimate 


m(k) 

(4.9) liminf — Tj < 1 — £ 

k—too k ^ J 
i= 1 

implies the existence of a subsequence Tj n satisfying 

(4.10) Ti + ■ ■ • + T in < ^4±i. 

£ 

Since ps,S',r' < (1 — p*) T / T *~ 1 for every positive integer t' and every S , S ' £ S, 
the probability that a sequence Ti contains at least one element T^+i = t' and, 
simultaneously, satisfies e{T\ +- T^) <T it+ \ = t' can be estimated as follows: 

P(T i+ 1 = r', Ti + • ■ ■ + Ti < Tj+i/e for at least one i) < —(1 — p*) T / ’ T * _1 . 

£ 

Therefore, 


P{Ti + 1 > t , Ti + • • ■ + Ti < Ti + i/e for at least one i) < -(1 — p*)^ * 

J = T 

As the right hand side of this expression tends to zero with t' —> oo and Tj n —> oo in 

(4.10), we conclude that relationship (4.9) holds with zero probability for each e > 0. 
This proves (4.8). 

Now, estimate (3.8) can be rewritten as 



Combining this formula with relations (4.2), (4.3) and (4.8), and taking into account 

that m{k) —> oo as k — » oo almost surely (because the Markov chain SV X 4 \-T k is, by 

assumption, well defined), we obtain (4.7). □ 

5. Numerical examples. In this section, we consider two examples of the input 
process and evaluate the convergence rate of the norm ||I?fc|| of the error for the 
iterative scheme (2.1), and that depending on the quantities T* and No as in (3.9) 
and (4.2), respectively. These are the two factors that affect the convergence rate 
guaranteed by Theorem 4.2. The following numerical computations were performed 
using the discrete dynamic Preisach (DDP) modeling algorithm which is suitable for 







14 


D. RACHINSKII AND M. RUDERMAN 


realization of the iterative scheme (2.1). The DDP algorithm has been introduced 
in [25] and applied for inverse hysteresis control in [26]. 

As the density function p(a, 0) of the Preisach model, we used the 5-finction 
p = 8{p — po) located at the point po = (a, (3) = (Ne/2,Ne/2). The convergence 
in this case is relatively slow, because the algorithm makes a nonzero correction to 
the approximating density p only when the switching region contains the point p o. 
Moreover, the correction is each time spread over multiple (at least (N — l)/2) nodes. 
The initial approximation of the scheme was the uniform density po(a, /?) = a with a 
small positive a. 

Two meshes with N l = 31 and IV 11 = 101 were considered (the odd number is 
conditioned by the numerical implementation of DDP algorithm). The number of 
nodes for these meshes, M 1 = 496 and M u = 5151, differs by one order of magnitude. 

5.1. Stochastic inputs. We will consider two types of the input Markov chain 
Xk- The first input denoted by is the Bernoulli (memoryless) process: 

P(x£=ie) = ^ l , i = 0 ,..., N. 

The second input process denoted by X® is a modification of the Bernoulli process 
where the sign of the input increment alternates at each time step: 

P(X» +1 =je\X*=ie) = \ J for odd fc; 

p i x k+i = J £ \ x k = * e ) = { i Q 1 ’ j > | for even k, 


where i, j = 0, ...,1V. Thus, a specific feature of this process is that it produces 
continuously oscillating input sequences. Note that both processes A and B have no 
knowledge of dynamics in the state space S of the Preisach model. 

We now show that the above two stochastic input processes satisfy the conditions 
of Theorem 4.2. 

Lemma 5.1. The Markov chain Yk = St h _ \-T k induced hy each of the processes 

A and B in the Preisach state space S is regular. 

Proof. A state S of the Preisach model will be identified with a polyline connecting 
the point (a,/3) = (0 ,Ne) with the line a = /? on the half-plane II. This polyline 
consists of vertical and horizontal links where the length of each link is a multiple of 
e. The relays are in state 1 below the polyline and in state 0 above the polyline [18]. 

The following finite sequence of input values will be called a standard sequence: 
(5.1) 


0, 

Xi, 









0, 

X2, 

Xi, 

0, 

X2, 






0, 

CO 

Xi, 

o, 

X3, 

X2, 

o, 

X3, 



0, 

Xi, 

Xu 

o, 

Xi, 

X2, 

o, 

Xi, 

X3, 0, 

X4 

0, 

X N , 

X 1 , 

o, 

Xl V, 

X2, 

0, 

Xi V, 

CO 

a 

0, 


0 , 


where Xi = is. This input sequence produces all the possible triangles (3.5). By 
definition, each of the processes A and B has a positive transition probability from 
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any state Xj to any state Xj. Therefore, for any given initial state, the input makes 
the standard sequence of steps and creates all the triangles with a positive probability. 
Hence, the process Yk = 5^+—+i>, in the Preisach state space is well defined (that is 
the times Ti are almost surely all finite). 

Now, it suffices to show that, given any pair of states S', S" £ S, the transition 

probability of the Markov chain 1*, = 5 th _t- T k from S' to S" is positive. Denote 

by Pq, ..., P' K , the corners of the polyline S' and by Pq, ..., P'^„ the corners of the 
polyline 5", where Pq = Pq' = (0,eN), a' K , = P' K ,, a'^„ = P'^,, and we use the 
notation P[ = (a', /3'), P" = (a", p"). For any state S" = Pq P" ■ ■ ■ P^„ £ S (except 
for the state 5" = {{a, P) : a = 0,/3 £ [0, Ne]} consisting of one vertical link), we 
define an input sequence associated with S'; this sequence is given by: 

• P 2 > a 3 > PI 11*5, • • •, oi'k"-i, P'k" if the link Pq'P" is vertical and the link P'k"-\P'k 
is horizontal; 

• P '^, ag, P'l, a",..., P'k"-\, a k" if both the links Pq P" and P'k"-\P'k" are 
vertical; 

• P",a' 2 , P 3 ,a'l,... ,oik"-i, Pk" if both the links Pq'P" and P'k"-\P'k" are 
horizontal; 

• p", ot'2, P 3 , a'l,... ,Pk"-i, if the link Pq'P" is horizontal and the link 
P'k"P'k" is vertical. 

The input sequence associated with the state S" induces the transition of the process 
5fc from the state consisting of one vertical link to the state S". 

Depending on a few features of the states S', S", we will define explicitly a finite 
sequence of input values, which induces the transition of the Markov chain Yk from 
the state Yk = S' to the state Yk +1 = 5"; the probability of such transition is positive, 
because the input sequence is finite. We call such an input sequence a transition input 
sequence. 

First, consider any state S' = PqP[ ■ ■ ■ P' K , which has a vertical link P'k'-iP'k" 
Consider the modified standard sequence (5.1) where we replace every instance of 0, 
except for the first and the last one, with the pair of entries 0,Xi: 

(5.2) 

0, xi, 


0, 

Xi, 

X2, 

Xi, 

o' 

Xi, 

X'2, 




cT . 

Xi, 

X.3, 

Xi, 

0 , 

Xi, 

X'3, 

Xi, 0, 

Xl, X3, 


0 ■ • 

Xi, 

X N , 

X 1 , 

0 , 

Xi, 

X N , 

x 2 , 0, 

N 

Xl, x'\ x 3 , 

■ , 0 , Xl, x N 


0. 

• The input sequence (5.2) induces the transition of the Markov chain Yk = 
S^+.-.+Tk from the state S' to the state S" = {( a,p ) : a = 0, p £ [0,iVe]} 
consisting of one vertical link. 

• If the polyline S" = Pq'P" ■ ■ ■ P^„ has more than one link and its link 
P’k"-\ ^k" is vertical, then a transition input sequence from S' to S" is the 
sequence (5.2), where we omit the successive four entries x\ ,P'k»-i,P'k>;0, 
which is followed by the input sequence associated with the state S". 

• If the polyline S" = Pq P” • • • P'k" has more than two links and its link 
P'k"- i P'k" i s horizontal, then a transition input sequence from S' to S" is the 
sequence (5.2), where we omit the successive four entries xi, a'^,, , a" K n_ x , 0, 
which is followed by the input sequence associated with the state S". 

• If 5" = PqP 1 = {(cx,P) : a £ [0,Ne],P = Ne} (i.e., S" consists of one 
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horizontal link), or if S" = PqP['P .has two links and the link P"P£ is 
horizontal and has the length of at least 2e , then a transition input sequence 
from S' to S” is the sequence (5.2), where we omit the successive three entries 
0, which is followed by the input sequence associated with the state 
S". 

• Finally, if S" = PqP”P% has two links and the link P”P!{ has the length e 
(i.e., P"P!J = {(a,/3) : a £ [0,e],/3 = e}), then a transition input sequence 
from S' to S" is the standard sequence (5.1), where we omit the first two 
entries, 0,iri, and add an extra entry X\ after the last 0. 

These cases define transitions with positive probability to all the destination states 
S" £ S from all the initial states S' = PqP[ ■ ■ ■ P'x’ with the horizontal link P'k'-iP'k 1 
. In each case, the transition input sequence creates all the triangles (3.5); the last 
triangle is created at the last step after which the destination state S" is achieved. 

The transition sequences from S' to S" are similar in the complementary case 
when the link P'k'-\P'k' °f the initial state S' is horizontal. In this case, a transition 
sequence consists of the values xjv,a,0 followed by the transition sequence defined 
above for the transition to the state S" (that is, we just add three entries at the 
beginning of the transition sequence). Here a is any input value satisfying a / 0, 
a ^ xn, a ^ □ 

5.2. Comparison of convergence rates. First we compare the convergence 
rates of ||£4|| for the processes X A and X B . Figure 2 presents simulation results shown 
by thick lines and exponential fits of the form ||£)-|| = cexp(— kS) (thin lines). For 
both input processes, the scheme (2.1) clearly demonstrates exponential convergence. 
The free parameters c and 5 of the fit are determined by using a standard nonlinear 



Fig. 2. Convergence of the mean squared error \\Ek\\ for the Markov processes X A and X B , and 
deterministic process X ( (thick lines). The simulation results are compared with the exponential 
fits ||-Efc|| = cexp(— k8) (thin lines). 

least-squares algorithm. Note that the exponential rate 5 is an empirical measure of 
the convergence rate A estimated from below in Theorems 3.1 and 4.2. According 
to Theorem 4.2, A > T**/(32iVoT*). We have estimated the time average T* of 
the times T and the time average T** of the inverse times 1/(1 + T*,) (cf. (4.2), 
(4.3)) for the simulations shown in Figure 2. We found that T** « 1/T*, hence 
A > n := 1/(32NqT^). In particular, the lower estimate fi of the convergence rate 
increases with the decrease of the average time T* required by an input process to 
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create all the triangles (3.5). Figure 2 shows that, in line with this observation, the 
actual empirical convergence rate <5 is also higher for the process A, which has smaller 
characteristic time T* than the process B. Table 1 compares the ratio 5 A /8 B of the 
empirical convergence rates for processes A and B with the ratio of the theoretical 
lower estimates ~ {T B ) 2 / (T A ) 2 of these rates guaranteed by Theorem 4.2. 

Table 1 

Evaluated, features of convergence rates for processes X A and X B . 


Input 

N 

rr tA 

J- * 

rr tB 

J- * 

5 a /5 b 

(T b ) 2 /(T a ) 2 

A a , X b 

31 

15778 

25934 

1.64 

2.70 


The convergence of the error of iterative scheme (2.1) for processes A' A and X B 
is further compared with that obtained using a deterministic input, denoted as X c , 
which constructs successively all the triangles (3.5). This is the same input sequence as 
used by the classical deterministic identification algorithm of Mayergoyz, but applied 
repeatedly as our algorithm is a recursive scheme. The input sequence is a concate¬ 
nation of inputs generating a mesh of first-order descending (FOD) curves [18,19]. 
The convergence rate 5 C for this input is higher than for stochastic inputs A a,b (see 
Figure 2) as expected, because the deterministic input has the shortest time T* (which 
is also deterministic in this case, T). = T]p for all k). 

Next, we compare the convergence rate for two different discretization meshes 
with TV 1 = 31 and IV 11 = 101 input states, respectively, and the corresponding meshes 
(3.1) on the Preisach operator half-plane II. Figure 3 presents simulation results 
obtained for the Bernoulli input process (type A). For each of the two meshes we 
observe the exponential convergence of the error. The convergence is much faster 
for the mesh with a smaller number of nodes, N 1 . This agrees with the fact that 
the theoretical lower bound of the convergence rate ft = T**/(32A i oT’*) is inverse 
proportional to Nq = 0(TV 4 ) (cf. (3.9)). Table 2 compares the ratio 5 1 /6 11 of the 
empirical convergence rates for the two meshes with the ratio of the theoretical lower 
estimates ft 1 /ft 11 ~ (TV 11 ) 4 /(TV 1 ) 4 of these rates guaranteed by Theorem 4.2. 





“ " " simul N 1 
expon N 1 


_ XT n 




0123456789 10 


Fig. 3. Convergence of the mean squared error of scheme (2.1) with the Bernoulli input (type 
A) for two meshes with N 1 = 31 and N u = 101 input nodes (thick lines) with exponential fits shown 
by thin lines. 
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Table 2 

Evaluated features of convergence rates for two meshes N 1 and A 11 . 


Input 

IV 1 

N n 

S'/S 11 

(N 11 f/{N 1 f 

X A 

31 

101 

16.40 

112.68 


6. Conclusions. The classical deterministic identification algorithm of Mayer- 
goyz recovers the density function of the Preisach model from input-output data that 
are obtained by testing the system with a series of deterministic test inputs which 
create a sufficiently dense mesh of first-order descending hysteresis curves. In this 
work, we have analyzed a recursive iterative identification algorithm that uses an 
input-output data stream generated by a random input process. This algorithm up¬ 
dates the values of the density function in the switching region at each time step. We 
have shown that the error converges to zero exponentially and obtained an explicit 
estimate of the convergence rate for a wide class of stochastic input processes. The 
rate of convergence guaranteed by this estimate depends on the target accuracy of 
the algorithm and on a certain stopping time T, which is a characteristic of the in¬ 
put process. Essentially, during the time interval of length T the input creates an 
input-output data set that contains sufficient information for recovering the density 
function. We have tested the recursive algorithm with two examples of stochastic 
input processes and the deterministic input sequence prescribed by the Mayergoyz 
algorithm and evaluated the rate of convergence numerically. We have found that the 
convergence rate decreases with increasing mean value of T and the number of nodes 
in the domain of the density function, which is in agreement with our theoretical 
estimate of the convergence rate. 
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